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Abstract 

To further improve the performance of Monte Carlo simulations of 
first-order phase transitions we propose to combine the multicanonical 
approach with multigrid techniques. We report tests of this proposi- 
tion for the d-dimensional $ 4 field theory in two different situations. 
First, we study quantum tunneling for d = 1 in the continuum limit, 
and second, we investigate first-order phase transitions for d = 2 in the 
infinite volume limit. Compared with standard multicanonical simu- 
lations we obtain improvement factors of several resp. of about one 
order of magnitude. 
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At first-order phase transitions standard Monte Carlo simulations in the 
canonical ensemble exhibit a supercritical slowing down JlJ. Here extremely 
large autocorrelation times are caused by strongly suppressed transitions be- 
tween coexisting phases which, on finite periodic lattices, can only proceed 
via mixed phase configurations containing two interfaces. Since the proba- 
bility of such configurations is suppressed by a factor exp(— 2aL d ~ 1 ), where 
a is the interface tension and L d ~ x the cross-section of the system, the au- 
tocorrelation times in the simulation grow exponentially with the size of the 
system, r oc exp(2erL d_1 ). 

A way to overcome this problem, known as umbrella 0] or multicanonical 
H sampling, is to simulate an auxiliary distribution in which the mixed phase 
configurations have the same weight as the pure phases and canonical expec- 
tations are computed by reweighting ||. Several tests for various models 
[0] have demonstrated that this method works well in practice and reduces 
supercritical slowing down to a power-like behavior with r oc V a = L da , 
where a ~ 1. While this is clearly an important step forward the remaining 
slowing down problem is still severe. In most cases it is even worse than for 
standard (e.g., Metropolis or heat-bath) Monte Carlo simulations of critical 
phenomena |J. 

For the latter applications several update algorithms have been developed 
which greatly reduce or even completely eliminate the critical slowing down 
problem |7]]. In addition to overrelaxation and cluster methods, an important 
class of such algorithms are multigrid techniques || ||]. Here the general 
strategy is to perform collective updates on different length scales by visiting 
various coarsened grids in a systematic, recursively defined way, generally 
known as V- or W-cycle |10 |. 

Because of their conceptual simplicity both the multicanonical reweight- 
ing approach and the multigrid update techniques are quite generally ap- 
plicable. The purpose of this note is to show that the two approaches can 
easily be combined and give a much better performance than each compo- 
nent alone. We report tests of this combination for the $ 4 lattice field theory 
with negative mass term in two conceptually different situations. We first 
consider the quantum mechanical tunneling problem in one dimension and 
study the performance of the new algorithm in the continuum limit. We then 
discuss field driven first-order phase transitions in the two-dimensional case 
and investigate the behavior of the multicanonical multigrid algorithm in the 
infinite volume limit. 
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For Potts models an interesting different approach was proposed only 



recently in Ref . |pJ|| . Here the idea is to combine a multicanonical demon 
algorithm with cluster update methods in a hybrid-like fashion. 

The basic idea of the multicanonical approach is to sample the mixed 
phase configurations with the same statistical weight as the configurations of 
the pure phases. At a field driven first-order phase transition this can always 
be achieved by a suitably chosen reweighting factor w _1 (m) = exp(— f(m)), 
where m = Y^i $i/V is the average field. In a temperature driven transition, 
m has simply to be replaced by the average energy. Starting from an initial 
guess based on experience or on some analytical approximation, a few iter- 
ations are usually sufficient to adjust this factor. Once it is fixed, canonical 
expectation values (0) Cim of any observable O can be computed from the 
basic reweighting formula 

(O) =^ (1) 

Wean ^ , W 

where (. . .) without subscripts denote expectation values in the multicanon- 
ical distribution. To update field values with a Metropolis algorithm in the 
multicanonical approach, we consider as usual local moves $j — > $j + A$j 
and compute the energy difference AE. The decision of whether such moves 
are accepted or not, however, is now based on the value of AE + f(m + 
A^/V)-f(m). 

The basic idea of multigrid techniques is to perform updates on dif- 
ferent length scales. Using the so-called linear interpolation scheme this 
amounts, in the equivalent unigrid viewpoint, to proposing moves for blocks 
of 1, 2 d , 4 d , . . . , V = L d = 2 nd adjacent variables in conjunction, with the 
sequence of length scales 2 k ,k = 0,...,n chosen in a specific, recursively 
defined order. Particular successful sequences are the so-called V-cycle with 
k — 0, 1, . . . , n — 1, n, n — 1, . . . , 1, and the W-cycle whose graphical repre- 
sentation looks like the letter W (for n=3, e.g., this is 0, 1, 2, 3, 2, 3, 2, 1, 2, 
3, 2, 3, 2, 1, 0, and for large n the W looks more and more like a "fractal"). 
In a canonical simulation the update at level k thus consists in considering 
a common move A$ for all 2 kd variables of one block, $j — > + A$, 
i G block. 

The modifications for a multicanonical multigrid simulation are quite triv- 
ial. Since at level k the proposed move would change the average field by 
2 kd AQ/V, the decision of acceptance is now simply to be based on the value 
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of AE + f(m + 2 kd A$/V) - f(m), with AE computed as in the canonical 
case. For didactic reasons we have emphasized here the conceptually sim- 
pler unigrid viewpoint. It should be stressed that in the recursive multigrid 
formulation, which can be implemented more efficiently (similar to the fast- 
Fourier transformation FFT), the multicanonical modification is precisely 
the same. 

We have tested the multicanonical multigrid algorithm for the scalar $ 4 
lattice field theory in d = 1 and d = 2 dimensions, defined by the partition 
function 



z = U 



d$i/A 



exp 



-^E(^(V^) 2 -^ 2 + ^)), (2) 



with A = \f2-ne and /U 2 ,g > 0. We always impose periodic boundary condi- 
tions. For d = 1 we keep Le = (3 fixed. Here the model describes the quantum 
statistics of a particle tunneling back and forth in a double-well potential in 
contact with a heat-bath at temperature T = 1/(3 |12| . At fixed (3 the limit 
L — ► oo corresponds to the continuum limit. For d > 2 we put e = 1. 
Here reflection symmetry is spontaneously broken for all fi 2 > fj^(g) > as 
L — > oo, which is now the infinite volume limit. Consequently, if a term 
h Si ®i is added to the energy, the system exhibits a line of first-order phase 
transitions driven by the field h. 

Even though in the one- dimensional case no spontaneous symmetry break- 
ing occurs, the numerical difficulties are quite similar. This is due to the 
fact that for small quartic coupling g tunneling events are strongly sup- 
pressed by a factor ~ exp(— const /g), which plays a similar role as the factor 
exp(— 2aL d ~ 1 ) at a first-order phase transition. The important difference is, 
of course, that the suppression factor stays roughly constant in the contin- 
uum limit, while at a first-order phase transition it rapidly decreases in the 
infinite volume limit. Nevertheless, for small values of g analogous slowing 
down problems in canonical simulations of the quantum problem are notori- 
ous, and a number of modified Monte Carlo schemes have been proposed in 
the past [|13|, [L4|]. None of these techniques, however, is general enough to be 
easily adapted to different potential shapes. 

To evaluate the performance of the multicanonical multigrid algorithm, 
we have recorded the time series for several observables and studied their 
autocorrelation times. In this note we shall concentrate on the average field 
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m = Yli^i/V) which reflects most directly the tunneling process. 

In previous investigations || emphasis was laid on the exponential au- 
tocorrelation time rffi of m, i.e., directly on the multicanonical dynamics. 
While this nicely illustrated the absence of exponential slowing down, it is 
not immediately clear how the remaining autocorrelations enter into the error 
estimates for canonical expectation values computed according to (p]) . To be 
precise, we are interested in the variance e 2 = cr^ = (0 2 ) — (0) 2 of the (weakly 

biased) estimator for (C) can , O = J2i m wim^Oi/ Y.\ m w{rrii) = WiOi/W i} if 
N m (multicanonical) measurements are performed. To facilitate a direct 
comparison with canonical simulations, we hence define for multicanonical 
simulations an effective autocorrelation time r cff by the standard error for- 
mula for N m correlated measurements, 

= vl n 2T^/N mj (3) 

where cr 2 an = (0?) can — (Ci)c an * s the variance of the canonical distribution of 
single measurements, which can be computed in a multicanonical simulation 
by using eq. ([!]). The squared error e 2 can be estimated either by block- 
ing or better by jack-knife blocking procedures, or by applying standard 



error propagation to the variance of O = WiOi/wi, which involves the (multi- 
canonical) variances and covariances oiwiOi and u>j, and the three associated 
autocorrelation times T m;m = T m , T wm - wm = T wm , and T wm - m = T m;U)m 16] . By 
symmetry, for O = m this simplifies to 



o (WiTnijWiTtli) 2,T wm 2 2t u , ?t . 



2 "'win 2 "'wm / a\ 



where (x;y) = (xy) - (x)(y) and r x . y = 1/2 + J2k( x o;yk}/(xo;yo) is the 
integrated autocorrelation time of multicanonical measurements. In this way 
properties of the multicanonical distribution (given by cr^J are disentangled 
from properties of the update algorithm (given by r wm ). Note that in r cff = 
( a muca/ a can) r «im, it is the autocorrelation time of w(rn)m that enters and not 
that of m, as previously investigated. 

Let us first discuss our results for the quantum mechanical case in d = 1, 
where we shall confine ourselves to the case /i 2 = 1.0, g = 0.04 and (3 = 10. 
This choice of parameters may perhaps be better characterized by the first 
few energy eigenvalues (obtained by a numerical integration of the associated 
Schrddinger equation), E = -0.913371, E x = -0.892348, E 2 = 0.029846, 



4 



and E 3 = 0.37813, or the probability ratio P^/ 'P max = P(0)/P max « 1.9 x 
10~ 3 , where P(m) is the probability distribution of the magnetization (or, in 
the present interpretation, of the average path). In a canonical simulation 
it would consequently be about 500 times harder to sample configurations 
with m = than configurations contributing to the peaks of P(m). To 
set up the multicanonical reweighting factor we started from a variational 
approximation jljj that works quite well for not too large /9-values and is 
known to provide locally a lower bound on P(m) ]Tj|. Alternatively, as the 
distribution depends only weakly on L, one could also use the distribution for 
small L (which is quite easy to generate by standard techniques) as input for 
the other simulations. Since the m-values vary continuously, we introduced 
bins of size Am = 0.02 to store the weight factor w(m). A single short run 
was usually sufficient to improve the initial guess such that the multicanonical 
distribution P'(m) had the desired flat shape, P'(m) ~ const, between the 
two peaks. 

In this way we performed multicanonical Metropolis and multigrid simu- 
lations for L = 4, 8, 16, 32, 64, 128, and 256, and using the multigrid update 
also for L = 512. In the multigrid case we investigated both the V-cycle 
(using n pre = n pos t = 1 pre- resp. post-sweeps || |9], [TOP and the W-cycle 
(using n pre = n post = 1 as well as n pre = 1, n post = 0). In the log- log plot 
of Fig. 1 we show our results for r eff of the multicanonical Metropolis and 
W-cycle (without post-sweeps) update algorithm, and for comparison also 
previous results |L9f for the canonical counterparts. We see that canonical 
and multicanonical simulations exhibit qualitatively the same behavior in the 
continuum limit L — > oo. For both distributions, the Metropolis update leads 
to a power-law growth r oc L z with z ~ 2, while for the W-cycle update the 
autocorrelation times stay roughly constant [pOfl . Here it is the overall scale 
which is reduced in the multicanonical simulation. For our choice of param- 
eters we obtain an improvement factor of about 60 (30) for the Metropolis 
(W-cycle) update. For smaller g, since this factor is essentially given by the 
inverse of the suppression factor exp(— const /g), we found the multicanoni- 
cal Metropolis update to be more favorable than the canonical W-cycle for 
reasonably large L. Eventually, however, there will always be a crossover at 
some L. For g = 0.04, if we compare the Metropolis update and the W-cycle 
in the multicanonical simulation we obtain an improvement factor of about 
50, 000 for L = 512. 

In the continuum limit the distribution P(m) and the ratio P m m/Pmax 
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stay roughly constant. This implies that the multicanonical approach can 
only give an improvement factor which is independent of L. In this case it is 
the update algorithm that plays the dominant role asymptotically for large 
L, while the multicanonical reweighting procedure sets the overall scale. The 
relative importance of the two approaches is just reversed at first-order phase 
transitions when the infinite volume limit is considered. 

As a test case we have studied the model @ with e = 1. For this model 
the line of second-order phase transitions separating the broken and unbro- 
ken phase in the /z 2 — g plane has recently been determined by Toral and 



Chakrabarti []2T|. Here we concentrate on the first-order phase transition 
between the two ordered phases at g = 0.25 and // 2 = 1.30, which is suffi- 
ciently far away from the critical point at /x 2 = 1.265(5) |2TJ to display the 
typical behavior already on quite small lattices. A sensitive measure of the 
strength of the transition is the interface tension a oa between the + and — 
phase, which turns out 0] to be a 00 = 0.03459(49) (for comparison, about 
the same value is found for the order- disorder interface tension in the two- 
dimensional g-state Potts model with q — 9, where a Q d = 0.03355. . . p2| ). 
We performed multicanonical simulations using the Metropolis update and 
the W-cycle without post-sweeps for lattices of size V = L 2 with L = 8, 16 
and 32. With the multigrid algorithm, due to the improved performance, we 
were also able to study lattices of size L = 64. Here each time series contains 
a total of 10 6 measurements taken every n e th sweep, after discarding 10 4 x n e 
sweeps for thermalization. The number of sweeps between measurements, 
n e , was adjusted in such a way that in each simulation the measurements of 
wm had an autocorrelation time of maximal 50, i.e., the length of each time 
series is at least 20, 000 r wm . 

A few of our results are collected in Table 1, where we give the integrated 
and exponential autocorrelation times of m and w(m)m as well as r cS ac- 
cording to eq. (|3|) for both update algorithms. We see that integrated and 
exponential autocorrelation times for m agree well with each other, show- 
ing that the corresponding autocorrelation function can be approximated by 
a single exponential. For wm we obtain values for that are consistent 
with those for m within error bars. The integrated autocorrelation times, 
however, are significantly lower, implying that the autocorrelation function 
is composed of many different modes. We also observe that the difference 
between r wm and r efr can be quite appreciable. From L = 8 to L = 64 the ra- 
tio r eS /r wm = ^muca/^can varies from about 1.9 to 4.6, reflecting the varying 
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probability distribution shapes with increasing L. 

By fitting r cff to a power law, r cff oc L z , we obtain for both update 
algorithms an exponent of z ~ 2.3, i.e., in this case it is thus the multigrid 
update that reduces the overall scale. The autocorrelation times of the W- 
cycle are reduced by a roughly constant factor of about 20 as compared with 
the Metropolis algorithm. Of course, for a fair comparison we should also 
take into account that a W-cycle requires more elementary operations than 
a Metropolis sweep ||. Such a work estimate, however, depends on many 
details of the implementation and it is hence difficult to give generally valid 
figures. With our implementations on a CRAY Y-MP we obtained a real 
time improvement factor of about 10. 



Table 1: Autocorrelation times for the multicanonical simulation using the 
standard Metropolis (M) or multigrid W-cycle (W) update algorithm. 







L = 8 


£=16 


L = 32 


L = 64 


T (0) 
m 


M 


212(12) 


668(23) 


3120(200) 






W 


11.30(32) 


37.2(2.0) 


148(11) 


746(62) 




M 


204.4(4.0) 


690(11) 


2984(63) 






W 


10.88(12) 


34.69(76) 


150.0(4.0) 


758(37) 


r (o) 

wm 


M 


209(12) 


655(31) 


2880(190) 






W 


11.34(33) 


36.9(2.0) 


146(13) 


600(120) 


T~wm 


M 


171.1(3.4) 


509.8(8.9) 


1840(40) 






W 


9.82(11) 


27.58(59) 


96.6(2.4) 


374(23) 


T eti' 


M 


322.7(6.1) 


1258(21) 


6050(120) 






W 


18.51(20) 


67.4(1.3) 


321.9(7.6) 


1724(86) 
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Figure Heading 



Fig. 1: Effective autocorrelation times t cS for the model (2) in d = 1 as 
a function of lattice size L with Le = (5 = 10, fi 2 = 1, g = 0.04 for 
different Monte Carlo algorithms. The canonical data are taken from 
Ref. P]. 
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